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We investigate the form of the viscous correction to the equilibrium distribution function in the 
context of a Cooper-Frye freeze out prescription for viscous hydrodynamic simulations of relativistic 
heavy ion collisions. The standard quadratic ansatz used by all groups for the case of shear viscosity 
is found to be disfavored by experimental data for Vi/{v2) 2 at the Relativistic Heavy Ion Collider, 
and is unlikely to be correct for the hadron resonance gas present at freeze out. Instead, data for 
V2(pt) along with V4,/{v2) 2 favor a momentum dependence between linear and quadratic. 



I. INTRODUCTION 



There has been much success in modeling the medium 
created in collisions at the Relativistic Heavy Ion Collider 
(RHIC) as a fluid. First the equations of ideal hydrody- 
namics were used with surprising success, indicating that 
the collision fireball appeared to behave like a nearly per- 
fect liquid [H-l|. More recently, this ideal hydrodynami- 
cal description has been generalized to viscous hydrody- 
namics in an effort to quantify dissipative corrections to 
perfect fluid behavior[6Hl0l|] 

The results of these simulations are, of course, sensi- 
tive to the initial conditions that are implemented for hy- 
drodynamic evolution, which are a topic of much study. 
However, the results can also be sensitive to the prescrip- 
tion used to describe the cessation of hydrodynamic be- 
havior once the system has become dilute ( "freeze out" ) . 
This freeze out prescription must describe the spectrum 
of particles that comes from a particular piece of fluid 
when it decouples. In principle this distribution is ob- 
tained from a kinetic description of the system at freeze 
out, which for a viscous fluid is of the form / = /o + Sf, 
where / is the equilibrium distribution and Sf is a small 
correction. 

The exact form of Sf, however, depends on the particle 
system in question. As pointed out in Ref. [Ill ], until 
recently (at least in the case of vanishing bulk viscosity), 
all groups have made the assumption that Sf oc p 2 fo- 
However, there is no particular reason to believe that 
this is the correct momentum dependence for the system 
of hadrons and resonances at freeze out in a relativistic 
heavy ion collision. 

In this article we use viscous hydrodynamic simulations 
to investigate the effect of using different forms for Sf on 
calculated observables in the absence of bulk viscosity 
and baryon chemical potential. We find that experimen- 
tal data, especially V4/ (W2) 2 , can be used to constrain Sf, 
in principle giving information about the dynamics of the 
hadron resonance gas present at freeze out. 



II. VISCOUS HYDRODYNAMICS 

If a system is near thermal equilibrium, it is useful to 
write its energy-momentum tensor as a sum of a locally 
isotropic (ideal fluid) part Tq V , and a viscous correction 
n A " y . If, following Landau, we define the local rest frame 
as the zero momentum frame T M ° = 0, we can write the 
result in terms of the four- velocity u M of this local rest 
frame with respect to some fixed lab frame: 

T „u = T »v +n ^ = £ u^u^-p A" v + tt^ + A^n . (1) 

Here = (g^ v — u^u v ) is a projector onto the space 
orthogonal to , while e and p are the energy density 
and pressure in the local rest frame. The viscous tensor 
W 1 " is symmetric and transverse to the fluid velocity: 
u ll H flL ' = 0. It is typically broken into a traceless shear 
tensor tt^ = (H^ — ^IY^g^ v ) and a term proportional 
to the trace, the bulk pressure n = |n^. To first order 
these terms are proportional to the shear viscosity rj and 
bulk viscosity £, respectively: 

ttM" = rjdbt-u") +0(d 2 ) (2) 

n = (d a u a + 0(d 2 ). (3) 

The equations of viscous hydrodynamics then follow 
from the local conservation of energy and momentum 

d^ v = Q . (4) 

Any other conserved quantity — such as baryon number — 
will have its own conservation equation and transport 
coefficients, but will be neglected here. Along with the 
equation of state of the particular system in question 
p = p(e) and the specified form of the viscous tensor (in 
principle including the temperature dependence of the 
transport coefficients), this uniquely specifies the evolu- 
tion of a viscous fluid. 

III. FREEZE OUT 

To model the breakdown of hydrodynamic behavior 
once a system has become cool and dilute, one must de- 
fine a freeze out process. This is done by taking a given 
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fluid cell (once it meets some predefined criteria such 
as reaching a certain temperature Tf) and converting it 
into a distribution of particles / given by kinetic theory. 
In the standard Cooper-Frye prescription, these particles 
undergo no further interaction, and the finally observed 
particle spectrum is given by integrating this distribution 
over the freeze out surface £ M : 
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where we have summed over the contribution from dif- 
ferent particle species i with degeneracy d;. Once thus 
obtained, it is often useful to express the angular depen- 
dence in terms of Fourier components, which for collisions 
of identical nuclei at midrapidity takes the form 
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For a system close to equilibrium, the distribution 
function can be written as the equilibrium Bose/Fermi 
distribution /o, plus a small viscous correction that van- 
ishes in the ideal hydrodynamic limit 



f(pn = fo(p-u)+sf( P ^. 



(7) 



The exact form oidf depends on the dynamics of the 
particles in question [ljj , so it will be useful to first derive 
the general properties that it must obey for any system. 
For simplicity, consider a single species with unit degen- 
eracy in the local rest frame of the fluid with 5f = 6f(p). 

For consistency, the system of particles should gener- 
ate the same energy-momentum tensor as the fluid just 
before freeze out: 



d 3 p 
(2n) 3 E 



which implies 
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The most general possible form, then, is: 
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Equation © for the spatial components L™ in the local 
rest frame gives a normalization condition for A and B 
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(11) 
(12) 



In the case of non- vanishing bulk pressure, Eq. @ gives 
yet another constraint. Recall that in the local rest frame 



we must also have u M II Ml/ 



n° M = 0. In particular, when 
we demand the Landau matching condition to be satisfied 
Il^uv = n 00 = 0, we find 
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which puts a nontrivial constraint on B(p) when II =/= 0. 

Note that taking B cx p 2 fo as in Ref. [12j fails to satisfy 
the Landau matching condition jl3[ — in fact, this is the 
case for any positive definite function. 

Investigating B(p) is left to future work. In the fol- 
lowing we will only concern ourselves with the case of 
vanishing bulk pressure, where B(p) plays no role. 

In contrast to B(p), there is no constraint on the func- 
tional form of A(p), which is related to the more standard 
notation x(p) by [11| 

A(p) = i/ (p)[l±/o(p)]x(p). (14) 

Until recently, all groups doing hydrodynamic simula- 
tions have assumed \ P 2 - Although this is correct for 
the most familiar solvable systems [ll(, there is no par- 
ticular reason to assume it to be correct at freeze out for 
a heavy ion collision, and we will see that that it is in 
fact disfavored by data. 



IV. MICROSCOPIC PERSPECTIVE 

In principle 5f is given by a kinetic description of the 
particles comprising the system at freeze out — here a col- 
lection of hadrons and resonances not far below the de- 
confinement transition temperature. In kinetic theory, 
the distribution function is determined by the Boltzmann 
equation 



(d t +u p -d x )/(p,x) = -C[/,p] 



(15) 



In equilibrium the collision term C[/, p] = 0, and the 
distribution function is 



/(P.x) = /o(p,x) = 
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By inserting / = fo + 5f into Eq. (fT5|) and keeping 
only terms to first order in df, one obtains the linearized 
Boltzmann equation. For a vanishing bulk viscosity, in 
the local rest frame of the fluid, this is I 111: 
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where C [6f, p] is the collision operator expanded to first 
order in Sf. Using Eqs. ©, (flU]) and (fbH) . and neglecting 
quantum statistics, one obtains after some algebra: 
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where <icr p k^p'k' is the differential cross section, P(9) 
(3 cos 2 9 — l)/2, and v re \ is the relative velocity: 
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Eq. ([IH]) is Eq. (B5) of Ref. [TTj] , with different notation 
to emphasize that one may also think in terms of cross 
sections rather than energy loss rates. 

The relevant cross sections here are hadronic cross sec- 
tions, since freeze out occurs in the hadronic phase. Let 
us estimate the typical energy of a hadronic collision at 
freeze out, which will indicate the relevant energy regime 
of these hadronic cross sections. In the above equation, 
p denotes the momentum of the hadron in the rest frame 
of the fluid. If the hadron momentum is parallel to the 
fluid momentum, p = u°pt — um t , where pt and mt are 
the transverse momentum and transverse mass in the lab 



frame, and u 
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The fluid velocity in realis- 



tic hydrodynamic simulations reaches values as large as 
u ~ 1. An observed pion with p t in the range 1-4 GeV 
comes from a part of the freeze out surface that has a fluid 
velocity that is parallel to its momentum and is near this 
maximum u. The momentum in the fluid frame p is sig- 
nificantly smaller — it is in the range 0.4-1.7 GeV; a col- 
lision with a thermal pion at T = 0.14 GeV then occurs 
at y/s in the range 0.4-1 GeV. Hadronic cross sections 
have a nontrivial dependence on y/s in this range. Ul- 
timately, a consistent description of heavy-ion collisions 
using viscous hydrodynamics should implement realistic 
hadronic cross sections, but this is beyond the scope of 
the present paper; the linearized Boltzmann equation is 
an integral equation that is nontrivial to solve even for 
relatively simple systems. 

However, one can gain insight on how xip) depends on 
the cross section by studying the limit where the momen- 
tum p is much larger than the temperature T. This limit 
is relevant to the study of anisotropic flow (v 2 and V4) 
at large transverse momentum p t (see Sec. IVI|) . In this 
limit, outgoing momenta p' and k' are essentially par- 
allel to p, so that P(9 PP >) ~ P(9 pk ') ~ 1, and Eq. ([T5|) 
simplifies to 
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If x(p) is a power law at large momentum p, xip) x P a 
with a > 0, one can neglect x(k) in the right-hand side. 
The energy of the collision in the center of mass is of 
order y/s ~ \/pP ■ By simple power counting, Eq. (|20j) 
then implies that the cross section varies with s like s 1 ~ a . 
The quadratic ansatz x{p) ^ P 2 solves the Boltzmann 
equation at large momentum only if the cross section 
decreases like 1/s. This is the case for, e.g, </> 4 theory [ill |. 
but not for hadrons in the relevant range of y/s, where 
hadronic cross sections rather increase with y/s. 



V. SETUP 

We would like to investigate the effect of using various 
forms for Sf at freeze out. To do this, we implement a 
2+1 dimensional conformal viscous hydrodynamic model 
with both Glauber- and CGC-type initial conditions, a 
constant value of rj/ s, a realistic equation of state, and a 
Cooper-Frye freeze out prescription including resonance 
feed down. All parameters are taken from the best fit 
results from Ref. [f|, to which we refer the reader for all 
relevant details. The only change here is to replace the 
quadratic ansatz for Sf (x = cp 2 ) with two other simple 
choices: a linear ansatz (x — cp), and a power law in 
between linear and quadratic (x = cp 15 ). It has been 
argued fill ] that the correct momentum dependence for 
Sf is expected to lie somewhere within this range. For 
each case, the same form of Sf is used for all particles at 
freeze out (using a Boltzmann equilibrium distribution, 
ignoring quantum statistics), with the overall tempera- 
ture dependent normalization c then fixed by the gener- 
alization of Eq. (JTTJ) : 
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Note that, in principle, even with the same momentum 
dependence, different particle species in a particular sys- 
tem could have different size contributions to the shear 
tensor (see, e.g., Sec. IV of Ref. [Hj]), but this possibility 
is not investigated here. 

All hadrons and resonances up to ~2 GeV are calcu- 
lated, with unstable resonances then allowed to decay 
after freeze out. 



VI. RESULTS 

Given the restriction of Eq. (fTTI) , any observable that is 
integrated over transverse momentum will be largely in- 
sensitive to the choice for Sf. We found this to be true in 
our calculations for, e.g., (pt) and multiplicity. However, 
the momentum dependence of differential quantities such 
as V2{pt) and v^ipt) can strongly depend on this choice. 
Figure [T] shows differential V2 and V4HV2) 2 for charged 
hadrons from our viscous hydrodynamic calculations as 
compared to experimental data from PHENIX at RHIC. 
As was found previously [ll|, Vi below p t ~ 1.5 GeV is 
largely unchanged for a range of reasonable choices for 
Sf. At larger p t , however, it can change the shape of the 
curve significantly. A stronger momentum dependence 
results in a larger viscous correction — and thus a larger 
deviation from ideal hydrodynamic behavior — at larger 
momentum. In contrast, Vi/(v2) 2 depends strongly on 
the choice for Sf even at smaller momentum. 

As can be seen in the plots for vi in Fig. [TJ the shape of 
the V2 curve depends on a combination of viscosity and 
the choice for Sf, while the normalization depends mostly 
on initial eccentricity and viscosity. Experimental data 
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FIG. 1. (Color online) Charged hadron V2 and Vi/(v2) 2 at impact parameter b — 8 f m for three choices of Sf: Linear (x oc p), 
Quadratic (x oc p 2 ), and in between (x oc p 1 ' 5 ), with ideal hydrodynamics results included for comparison. The viscous Glauber 
calculation is for rj/s = 0.08 and Ti = 0.333 while the CGC calculation used ri/s = 0.16 and Ti = 0.299 [f|. Experimental data 
for 1)4/(1)2) 2 scaled down by a factor of 1.37 to account for fluctuations [3, 0- Data at 30-35% centrality from PHENIX with 
statistical error bars only [la ]. 



appear to favor a comparatively large initial eccentricity, 
such as can accommodate a viscosity of order ij/s ~ 0.16, 
while the favored Sf depends on this viscosity value. A 
linear ansatz seems unlikely, while a stronger momentum 
dependence works well. 

The curve for ^/(t^) 2 , on the other hand, favors a 
weak momentum dependence. The quadratic ansatz in 
particular is quite difficult to reconcile with the very flat 
experimental data — at least for 77/s > 0.08 and/or for 
more peripheral collisions. A Sf with a stronger momen- 
tum dependence also results in a strong impact parame- 
ter dependence, which is not seen in data (see Ref. 15!]), 
while a linear ansatz results in a dependence on impact 
parameter that is even smaller than is seen from ideal 
hydrodynamics. 

Note that for any choice of Sf, the magnitude of 
va/ {V2) 2 is somewhat smaller than the experimental data, 
even after correcting for the estimated effect of eccentric- 
ity fluctuations and especially for CGC-type initial con- 
ditions. This could be due to the presence of additional 



fluctuations or differences in the initial energy density 
profile, and this is not yet fully understood. The shape 
of the curve as a function of pt, however, is not sensitive 
to these aspects of the calculation and is largely deter- 
mined by the viscous correction (in conjunction with the 
freeze out temperature) [l5j]. In addition, both STAR 
and PHENIX find a flat pt dependence, though at dif- 
ferent values. It therefore appears that systematic errors 
in the experimental data may affect the magnitude of 
V4,/(v2) 2 i but arc unlikely to strongly affect the shape. 
Thus, it is mostly the shape of the curve that matters 
here, and this is what provides the constraint on the vis- 
cous correction, clearly showing a preference for a mo- 
mentum dependence that is weaker than quadratic. 

Thus, data as well as theoretical considerations point 
to something between the linear and quadratic ansatz. 

One should refrain from trying to extract a precise 
value from the finer details of Fig. [TJ however, since cer- 
tain details could change with an improved treatment. 
For example, our model uses a simple Cooper-Frye freeze 
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out prescription with simultaneous chemical and kinetic 
freeze out, and the initial conditions are still quite un- 
certain (hence the use of two different types of initial 
conditions), and so not all the physics is perfectly cap- 
tured. Even within our model, an interplay with freeze 
out temperature could change which exact Sf is most fa- 
vored by W4/(f2) 2 data 15]. In addition, hydrodynamic 
results in general have uncertain applicability at large 
Pt- As p t increases, there is an increased probability that 
other physics is coming into play and one should have 
less confidence in hydrodynamic calculations. 



Nevertheless, it is clear that the uncertainty involved 
in the correct form for Sf has potentially measurable ef- 
fects. Further improvements in hydrodynamic calcula- 
tions could therefore potentially give valuable informa- 
tion about the microscopic dynamics at freeze out by 
comparing to data. Conversely, theoretical progress from 
the microscopic side could give information about where 
the validity of the hydrodynamic description starts to 
break down. 



VII. CONCLUSION 

It has been thus far unclear what is the correct form 
for the viscous correction to the equilibrium distribution 
function used at freeze out in viscous hydrodynamic cal- 
culations of heavy ion collisions. We investigated the 
effect of various reasonable choices for Sf in addition to 
the standard quadratic ansatz on the resulting differen- 
tial elliptic flow V2(pt) as well as the higher-order har- 
monic Vi{pt) and compared to experimental data. The 
quadratic ansatz is strongly disfavored by RHIC data, 
with the most likely form lying somewhere between lin- 
ear and quadratic. Specifically, data for va/(v2) 2 favor a 
weak momentum dependence in Sf both from the shape 
of the curve as a function of pt , as well as the centrality 
dependence. Data for V2(pt) favor a stronger-than- linear 
momentum dependence in Sf, but the precise best-fit 
choice is viscosity dependent. 
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